% This function runs the second step of the SR approach
function [theta2Step2Opt,AVarTheta2Step2] = step2SR(AVarTheta1Step1,outputStep1,bandwidthT,numBootStep2)

% option parameters
dispOn        = 1;
allowDeSaling = 1;           % allow for descaling of the var_u to ensure omegeHat is positive definite
seedNum       = 1;           % seed for the simulated paths in the bootstrap

[theta2Step2Unadj,deScaling,~,~,hxHat] = Var1_GMM_closedForm(AVarTheta1Step1.VarFactorsRobust,AVarTheta1Step1.autoCov10Factors,...
    AVarTheta1Step1.autoCov01Factors,outputStep1.factors,allowDeSaling);
if deScaling ~= 1 
    disp(['The measurement error was descalted by ', num2str(deScaling), ' in the estimated P dynamics'])
end

% Note: The market prices of risk are given by
%h0  = h0Q + f0;     %note MPR is scaled by sigma
%hx  = hxQ + fx;

% Adopt bias-adjustment
if numBootStep2 > 0
    biasCorrTheta2On = 1;
else
    biasCorrTheta2On = 0;
end
if biasCorrTheta2On > 0
    theta2Start      = theta2Step2Unadj;
    [theta2Step2Adj,errorMes]  = biasCorrectTheta2(theta2Start,AVarTheta1Step1.VarFactorsRobust,...
            AVarTheta1Step1.autoCov10Factors,AVarTheta1Step1.autoCov01Factors,...
            outputStep1.factors,bandwidthT,allowDeSaling,numBootStep2,seedNum,outputStep1.numBootStep1And3);
    
    % Reporting results
    if errorMes == 0
        % The bootstrap was succesful
        theta2Step2Opt            = theta2Step2Adj.theta2Adj;
        % The standard errors for theta2 using asymptotical inference
        AVarTheta2Step2 = getAVarTheta2(theta2Step2Opt,AVarTheta1Step1.VarFactorsRobust,...
            AVarTheta1Step1.autoCov10Factors,AVarTheta1Step1.autoCov01Factors,outputStep1.factors,bandwidthT);
        AVarTheta2Step2.names                   = fieldnames(theta2Step2Opt);
        AVarTheta2Step2.Draws                   = theta2Step2Adj.theta2BootDraws(:,1:outputStep1.numBootStep1And3);
        AVarTheta2Step2.VarBoot                 = theta2Step2Adj.VarBoot;
        AVarTheta2Step2.AVarBootDraws           = theta2Step2Adj.AVarTheta2BootDraws;
        AVarTheta2Step2.theta2_pTile            = theta2Step2Adj.theta2_pTile;
        AVarTheta2Step2.tTestTheta2_pTile       = theta2Step2Adj.tTestTheta2_pTile;
        AVarTheta2Step2.AbstTestTheta2_pTile    = theta2Step2Adj.AbstTestTheta2_pTile ;
    else
        % The bootstrap was NOT succesful and we simply use the OLS estimate
        theta2Step2Opt            = theta2Step2Unadj;

        % The standard errors for theta2 using asymptotical inference
        AVarTheta2Step2 = getAVarTheta2(theta2Step2Opt,AVarTheta1Step1.VarFactorsRobust,...
            AVarTheta1Step1.autoCov10Factors,AVarTheta1Step1.autoCov01Factors,outputStep1.factors,bandwidthT);
    end
else
	% We use the OLS estimate but ensure that the VAR model is stationary
    eigValues           = eig(hxHat);
    realEigValues       = real(eigValues);
    imagEigValues       = imag(eigValues);
    modulus             = sqrt(realEigValues.^2 + imagEigValues.^2);
    if any(modulus>=1)
        [h0HatScaled,hxHatScaled,sigmaHatScaled,scalingOpt] = runInduceStatVAR(hxHat,outputStep1.factors,AVarTheta1Step1.VarFactorsRobust,...
            AVarTheta1Step1.autoCov10Factors,AVarTheta1Step1.autoCov01Factors);
        theta2Step2Opt = setTheta2(h0HatScaled,hxHatScaled,sigmaHatScaled);
        if dispOn == 1
            disp(['The VAR estimates are scaled by ', num2str(scalingOpt),' to induce stationarity']);
            disp(['Unscaled eigenvalues = ',num2str(modulus')]);
            disp(['Scaled eigenvalues   = ',num2str(eig(hxHatScaled)')]);
        end
    else
        % No scaling needed of the VAR model
        theta2Step2Opt = theta2Step2Unadj;
    end
    % The standard errors for theta2 using asymptotical inference
    AVarTheta2Step2 = getAVarTheta2(theta2Step2Opt,AVarTheta1Step1.VarFactorsRobust,...
        AVarTheta1Step1.autoCov10Factors,AVarTheta1Step1.autoCov01Factors,outputStep1.factors,bandwidthT);
end


% Additional output
AVarTheta2Step2.deScaling = deScaling;

end

